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ABSTRACT 

We discuss two methods for reducing the variance of 
estimators of parameters of limiting distributions of stable 
stochastic processes in simulations. The methods are discussed 
in the context of the simple GI/G/1 queue. Of the two methods 
one, which we call an internal control variable, gives a vari- 
ance reduction which is roughly uniform over values of the 
parameters of the process and, in particular, works well for 
values of the traffic intensity, close to 1. 
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VARIANCE REDUCTION FOR REGENERATIVE SIMULATIONS, I: 

INTERNAL CONTROL AND STRATIFIED SAMPLING FOR QUEUES 

by 

Donald L. Iglehart and Peter A. W. Lewis 

1 . Introduction and Summary 

Simulators are frequently faced with the task of estimating 
a parameter associated with the limiting distribution of a stochas- 
tic process which is being simulated. A methodology based on 
regenerative processes for obtaining point estimates and confidence 
intervals for such parameters has recently been developed in Crane 
and Iglehart (1974a, b) , (1975a, b) and Iglehart (1975), (1976a, b) . 

In this paper we shall indicate two techniques, internal control 
variables and internal stratified sampling, which might be used in 
conjunction with the regenerative method for obtaining additional 
variance reduction for the estimates. To illustrate these tech- 
niques we shall restrict our attention in this paper to the GI/G/1 
queue. Other applications of these ideas will be dealt with in 
future publications, as well as uses of the regenerative property 
which may possibly be more suited for obtaining variance reduction 
for estimates in stable stochastic processes. 

In the GI/G/1 queue we assume the zeroth customer arrives 

at time tp = 0 , finds a free server, and experiences a service 

time Vq . The n — customer arrives at time t and experiences 

a service time v . Let the interarrival times t - t = u , 

n n n-1 n 

n>l. Assume the two sequences {v : n>0} and {u : n>l} each 

n n 



consist of independent, identically distributed (i.i.d.) random 

variables (r.v.'s) and are themselves independent. Let E{v }= 

n 

y~\ E{u n }=A - '*', and p = A/y, where 0<A,y<°°. Thus y has 

the interpretation of the mean service rate and A has the inter- 
pretation of the mean interarrival rate. The parameter p is 
called the traffic intensity and is the natural measure of conges- 
tion for this system. We shall assume that p< 1, a necessary 
and sufficient condition for the system to be stable. 

While many characteristics of interest can be estimated 

using the regenerative method, we shall restrict our attention to 

th 

the waiting time of the n — customer, W^ (time from arrival to 
commencement of service) . To obtain a representation for the pro- 
cess {W : n > 0} let X=v ,-u and set S A = 0, S =X. + ...+X , 
n n n-1 n 0 n 1 n 

n>l. The following well-known recursive relationship exists for 

the W ' s : 
n 



W_ = 0 



W 



n+1 



= [W +X , , ] 
n n+1 



n > 0 . 



By induction, we also have 



W = max{S -S, : k = 0 , 1 , . . . ,n } , n>0. 

n n k. 



Using the strong Markov property of the process {S^: n>0} 

it can be shown that there exists a sequence of integer-valued r.v.'s 

(Bj,: k > 0) such that Sq = 0, B^, < B^-^, and W^ = 0 with proba- 

p k 

bility one. In other words, the customers numbered B^ are those 
lucky fellows who arrive to find a free server and experience no wait- 
ing in the queue. The fact that there exists an infinite number of 
such customers in the GI/G/1 queue is a direct consequence of the 
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assumption that p<l and the strong law of large numbers. If 
we let a k = ' k-l/ then represents the number of 

customers served in the k— busy period (b.p.) and they are 
numbered . . . ,3^-1) . 

Next define the sequence {Y^: by 



V 1 



Y k = 



j=B 



k-1 



V 



th 



is stable for p < 1 we have W n => W as n-*- 00 , where 



the sum of the waiting times in the k — b.p. Since the queue 

» denotes 
n 

weak convergence. 

Our goal is to estimate E {W} by simulation . 

It is known that E{W} < 00 if and only if E{(X*) 2 } < °°. We 

shall, for simplicity, assume that E{v 2 } < 00 which will guarantee 

that E{w} <°°. The regenerative simulation method is based on the 

analytic results that the sequence {(Y^,a^): k - 1} i- s independent and 

identically distributed and E{W} = E{Y^}/E{a^} , the ratio of two 

means. This suggests using the ratio of estimates of E(Y ) and 

1 n 1 

E(a,) to estimate E (W) . Thus, if we let Y (n) =— £ Y, and 

1 1 n n k=l K 

a(n) =— £ a, , where n now denotes the number of cycles observed, 

n k=l k 

then a ratio estimate of E(W) , for example, is 



W(n) = Y (n) /a (n) . 



( 1 . 1 ) 



(In the sequel we drop the dependence on n unless necessary.) 

Now let Z^ = Yk - E{W}ak / k>l, and note that the sequence 
{Z^.: k>l} is i.i.d. and E{Z^} = 0. We assume the variance of Z^. , 
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E{z£} = o 2 = var{Y k >-2 cov{ Y R , a k >E {W> + var{a k )E 2 {W} (1.2) 

is finite and positive. Then one can easily show that as n -*■ 00 



n^/2 j-y y a ^ ^ -e{w} ] 
a/E{a j 



N (0 , 1) , 



(1.3) 



where N(0,1) is a mean zero, variance one normal random variable. 
This result yields a confidence interval for E{w} provided we 
can estimate a/E{a^}. A variety of estimates have been studied 
and are reported on in Iglehart (1975a) . Here we just mention two. 
The so-called classical estimate for a/ E{a^} is given by 

S 1 = [S lx - 2S 12 (Y/a) + S 22 (Y/a) 2 ] 1/2 /a. 



where S, n is the sample variance of the Y. ' s , S„ n of the a ' s , 
11 c k 22 k 

and S ^2 the sample covariance of the Y k ' s and a k 's. 

The second estimate of a/E{a 1 > is the jackknife estimate 
which is defined to be 

S= [ I (9 • -9) 2 / (n-1) ] 1/2 , 
z i=l 1 

- - i n 

where 0. = n(Y/a) - (n-1) ( £ Y ./ £ a. ) , 1 < i < n, and 0 = - Y 0 . 

j^i 3 j/i 1 * n J-=l 

is called the jackknifed estimator of E (W) . Both and S 2 

are strongly consistent. 

The jackknifing technique is known to work particularly 
well as a confidence interval estimate for ratios; for a large 
number of cycles n the computational problem is severe, but in 
that case the technique using (1.3) and works well. For 

details see Miller (1974) and Iglehart (1975). 
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I 



The problem addressed in this paper is how to apply variance 
reduction techniques with the ratio estimator W(n). In almost 
all practical situations, where in particular one might want to 
compare the mean waiting times of two different queueing systems 
(Iglehart (1976) ) , there is a premium on achieving the minimum 
possible variance for the estimation in the given computing time 
(number of cycles allowed) . Variance reduction techniques for 
simulations are discussed in Kleijnen (1974) and Gaver and Thompson 
(1973) , but there are difficulties in applying these to ratio 
estimates and regenerative systems. In particular variance reduc- 
tion via the usual control variable techniques is difficult. A 
variation of this technique, which we call internal control variables 
and which is generally useful for ratio estimates, is introduced 
and shown to give considerable variance reduction for the point esti- 
mate W(n). Another technique, internal stratified sampling, is also 
explored. It is a natural technique to use but appears to be diffi- 
cult to use with ratios. Moreover, it becomes less and less 
effective as p-»-l, while the internal control technique holds up 
well for p close to 1. In fact the internal control variables 
described here for an M/M/1 queue give a variance reduction which 
is fairly uniform for all values of y and ratios p<l. Better 
results can be obtained for particular cases of the parameter p. 

It will also be apparent after the development of the vari- 
ance reduction techniques in the next sections that the two tech- 
niques for confidence interval estimates discussed above apply to 
the estimates after variance reduction has been applied. 
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2. Internal Control 



Two main methods for variance reduction, antithetic variates 
and control variables (Gaver and Thompson (1973)) have been put 
forward for use with queues in the non-regenerative situation. Of 
these, antithetics has very limited utility beyond the simple M/M/1 
situation in which it is patently clear how to generate two reali- 
zations of samples which give negatively correlated estimates. 

That scheme for generating negatively correlated estimates does 
not work with the regenerative method because the regenerative 
blocks in the original realization of the queue and the antithetic 
realization get out of synchronization. No alternative way has 
been found to utilize antithetic variates with the regenerative 
technique and we feel, along with many computer scientists, that 
the technique has limited validity in systems simulation. 

The technique of using a control variable and, in particular, 
a regression-adjusted control variable (Gaver and Thompson (1973) 
p. 587) is much more broadly applicable in systems simulation, 
although it is again difficult to adapt to the regenerative situa- 
tion. Briefly, say we are simulating an M/G/l queue to estimate 
E(W) with the non - regenerative te chnique of averaging the first 
m waiting times to obtain an estimate of the unknown E(W), 







1 m 
- I 

m ^ 



W . . 

m j=l 3 



( 2 . 1 ) 



The same random numbers used to generate the m non-exponentially 
distributed service times are used to generate m exponential 
service times for simulating an M/M/1 queue whose input stream 
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is identical to that of the simulated M/G/l queue. One would 
then expect that if the G-distribution is not too different from 
the exponential distribution, successive waiting times, say , 
in the M/M/1 queue realization would be close to (positively 
correlated with) the corresponding waiting times in the 

M/G/l queue. Consequently, the average of the W('s, say W^, 
will be positively correlated with and one can form a new 

estimate 



W = W + B(W'-E(W')) • 
mm mm 



( 2 . 2 ) 



Now E (W^) is close to E{W' } = p/y (1-p) for m large, so 

E (W ) ~E(W ) and the variance of the new estimate will be a minimum 
m m 

if 



$ = -cov(W ,W' ) /var (W' ) ; 

mm m 



( 2 . 3 ) 



in fact 



o& & 

W s.d. (W ) 

m m 



ry rvy 

W s.d. (W ) 
m m 



- <l-r*) 1/2 . 



( 2 . 4 ) 



where 



^ ~ cov(W , W ) 

r = corr (W ,W ) = m m 

m m 



a~ 



W 
m m 



( 2 . 5 ) 



There are a number of important points to be noted about 
this technique: 

(i) It allows one to use known analytical results (such as 
the expected value of the limiting waiting time in an 
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M/M/1 queue) to reduce the variance of a simulation. 

Such use of analytical results is a basic principle of 
simulation. 

(ii) It can be extended to non-linear controls or to multiple 
control by more than one control variable. 

(iii) The great art in the technique lies in finding a control 
whose mean value is known and which is highly correlated 
with the estimator which is being controlled. Thus (2.4) 
says that |r| must be close to 0.9 to reduce the stand- 
ard deviation a~ to one half of o~ and this 

w w 

m m 

generally is equivalent to reducing the required sample 
size for a given precision by a quarter. Controls which 
are that highly correlated with the estimate are not 
easy to find. 

(iv) It is in general too much to ask that the correlation and 
variances in (2.3), (2.4), and (2.5) be known analyti- 

cally. They, therefore, must be estimated from the 
simulation data and this will reduce the variance reduc- 
tion which is attained. 

Now using this method to control the regenerative estimate 
given in the introduction , 

W(n) = Y(n)/S(n) , 

where n refers to a fixed number of cycles in the queue, one runs 
into the same problem of synchronization as with the antithetic 
case, namely n cycles in the M/G/l queue may take a very different 
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number of waiting times to achieve than are needed for n cycles 
in the M/M/1 queue. Moreover, the correlation between Y£ and 
Y^ will be weak and made even weaker by the use of the ratio 
estimate. This problem becomes rapidly apparent in simulation 
studies of the technique. 

The following technique which we have called internal 
(within block) control has been developed to overcome this. It 
is, in fact, a special case of the very broad technique called 
concomitant variables in Gaver and Thompson (1973) , p. 588 and 
will be illustrated only for point estimation of E(W) and E(a). 

Its extension to other quantities discussed by Crane and Iglehart 
(1974a) is immediate in principle, although the control quantities 
discussed below may be different. 

2 . 1 Internal Control: Basic Ideas 

The idea of an internal control variable is simple. In the 
estimate W(n) , the averages Y (n) and a(n) contain n random 
variables Y^ and which are each functions only of the 

'fch 

interarrival and service times occurring in the k — cycle (or b.p.) 

and are independent of the other interarrival times and service times. 

Thus, it is natural to use some function of these 2a^ random 

variables to control each Y. and a, . 

k k 

The naive application of this idea is that if we can reduce 
the variance of both the numerator and denominator Y(n) and a(n) 
we will reduce the variance of W(n) , but we will show that the 
situation is more complex than this. We will denote a function of 
the random variables in the k— cycle by C(k) but will generally 
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drop the index. In general we also use C T (k) to denote a control 
for the numerator (top) in the ratio estimator and C D (k) to 
denote a control for the denominator (bottom) . 

Typically, C(k), or simply C, could be the difference 



between the service time v 



6 , 



and the time to arrival of the next 



customer from the arrival of the Bj,th customer, namely 
This difference has, in a GI/G/1 queue, a known mean 



k+1 



y ^ - X and large positive values of this function C(k) corre- 
spond to large values of and a k , and vice versa. We return 

to specific control variables and their computation in the next 
sub-section . 

Note that one can control either the top or bottom of the 

ratio estimator, or both; to fix ideas assume we control the top 

and have, in general, an internally controlled estimate 

n 



W CT (n) = 



n ^ Y k + 6 T* C T E(C T ))} 



( 2 . 6 ) 



a 



where 6 T / as in the usual control estimation technique, is fixed 
so as to minimize var(W CT (n)). In practice it is usually esti- 
mated from the simulation data. 

Now the quantity a 2 / E 2 {a}, where cr 2 =E{z£} is given at 
(1.2) , is just the leading term in the asymptotic expansion for the 
variance of the ratio estimator W(n), and carries over to the 
more complicated situation (2.6) to give (asymptotically as n -*• °°) 



n var(W CT (n)) + 



2 

°T = {fro-} (var(Y'+6 T C^,-a') 

2 

= {Hff} var { (Y ' “Ot ' ) + B t C' t }, (2.7) 
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where Y'=Y/E{Y}, a'=a/E{a}, and C^,= C t /E{Y}. Fora 

derivation of this result see Cram4r (1946), p. 354, eq. (27.7.3). 

It follows from (2.7), as before, that to minimize the 
asymptotic variance of W CT (n) we must take 

cov(Y'-a' ,C' ) cov (Y ,C_) „, v . cov(a,C ) 

_ R = £_ - £__£d£i I— to q\ 

P T var(C,p var(C T ) E (a) var(C T ) 

But most importantly we notice that C^, must be highly correlated 
with the difference Y' - a'. To achieve this is much more difficult 
than finding a quantity which is highly correlated with either of 
Y or a, simply because Y and a are highly correlated and 
increase together. In particular if a=l, which has a high 
probability if p, the traffic intensity is small, then Y = 0. 

Note, too, that E (W) =E(Y)/E(a), the quantity we are trying 
to estimate in the simulation, appears in (2.8) for the top control 
regression coefficient. Also similar equations to (2.7) and (2.8) 
pertain to the case where the bottom of the ratio is controlled 
(in this case a) , and simultaneous equations can be derived for 
8™ and 8 if both top and bottom are controlled. The additional 
complexity in estimating 8 m and &_, does not seem to be justi- 
fied by simulation results (discussed later) which also show that 
if only one control is used there seems little to choose in terms 
of reduction achieved between putting it on the top of the bottom. 

In both cases the control must be highly correlated with the 
difference between Y and a. 
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The above results are generally applicable for any regenerative 
process in which ratio estimates are used. The choice of C is, 
of course, the art in the design of a simulation with variance 
reduction and is considered for the GI/G/1 and specifically the 
M/M/1 queue in the next section. In all cases this choice is 
limited by one ' s ability to compute analytically E (C) . 

2 . 2 Internal Control: Design Considerations 

We discuss here the design of internal controls for the 
M/M/1 queue, the ideas being applicable to the GI/G/1 case with 
the proviso that the computations might be considerably more diffi- 
cult. Simple computations of E(C) are given here and more diffi- 
cult ones in the Appendix; we do not distinguish between bottom 
and top controls, since both must be correlated with the difference 
Y-a, and we drop consideration of cycle number, since all variables 
are within the cycles which have identical structure. 

Again, we are considering estimation of E (W) , but of the 
many possible controls, those listed below would probably work as 
well with other functions of W, e.g. percentiles. The controls 
are listed roughly in order of complexity of computation and of 
supposed correlation with Y-a. This can usually only be guessed 
at and generally the more elaborate controls which might have 
greater correlation with Y-a are more difficult computationally. 

Superscripts on C are labels to differentiate the controls. 

We have discussed above the difference = Vq - u^ , whose moments 
are simple to compute. Then we have 
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0 



if 



a 



1 



(2.9) 



C 



( 1 ) 





= X ± if a > 2 . (2.10) 

It is easy to compute E(C^) for M/G/l queues and 
possible for the GI/G/1 queue. Thus, for the M/M/1 case we 
have, using the Markov property of the exponential distribution. 



P{v 0 <u 1 > = P{a=l) = | [1-F u (y) ]dF v (y) 



e _Xy (ye yy )dy = (2.11) 



1+P 



which goes from 1 to 1/2 as p = A/y goes from 0 to 1; furthermore. 



P f a > 2 } = 1 — i ^ 

Fl0t ~ Z1 1 1+p 1+p * 



( 2 . 12 ) 



Now given that x i = v o ~ u l 9 reater than zero, the 

excess v^ - u^ is distributed as an exponential random variable 
with parameter y. Therefore 

E(x i> =E(C (1) ) = o x ^ + t£p k - • < 2 - 13) 

The variance of can also be computed. 

One would generally like to obtain more correlation of C 
and Y - a when Y is large, and one feels this can be done by 
bringing in the second waiting time. Thus we have as control 
candidates 
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(2) 



(3) 



= 0 , 



x t + x 2' 



<X* + X 2 ) + 



w 2 ; 



if a = 1 
if a > 2 ; 



(2.14) 



(2.15) 





c <4, . 


_i_ j. 




if 


a 




x t 


+ (x 1 +x 2 ) , 




if 


a 


The 


(2) 

control C 


is W Q = 0 


if 


a = 


If 


a = 2 , and 


it is Wq + W^ 


+ X+ if a > 


3. 


It 


is 



(2.16) 



if 



the effect of the first two waiting times without the additional 
computational complexities involved in computing the expectations 



of 



(3) 



and 



(4) 



Simulation results show that C 

( 2 ) 



(3) 



(4) 



little more control than C 
we have 



and C ' give very 
for which, in the M/M/1 case. 



E(C (2) ) = 0 + E (X++X+ ,a>2) 





= 


p 

1+p 


{e(x++x 2 I 


x+> o> 








= 


p 

1+p 


& +E(X 2 1 


X 1 > °>} 








= 


p 

1+p 


& +E(X 2>} 


_ p a p n 

l+pjy 1+p y 


\ 


(2.17) 


using (2.12) 


and the 


fact 


that X 2 


is independent 


of X+. 


We use 


the notation 


e{x,a) 


for 


e{xi a >. 


where X is a 


random 


variable , 



A an event, and 1^ the indicator function of A. 

Similar computations go through for and C^; from 
these we will need later the following illustrative result (M/M/1 
case) . 
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P{a=2 } = P{X^>0,X^+v 1 <u 2 } 

= P{X^>0,v 0 -u 1 +v 1 <u 2 } 

= — {p ( u > Y } = — H s P 

l+p 1 ( 2 1 1+P (y+A) z (l+p) 3 ' 



(2.18) 



where Y, the sum of two exponential (y) random variables, has 
a Gamma (y, 2) distribution. 

From this we also get that 



P { a> 3 } 



1 P 

l+p (l+p) 3 



P 2 (2+p) 
(1+P) 3 * 



(2.19) 



Now, none of these four controls is specifically designed 
to be correlated with the difference Y - a. As a result they work 
well for y and A such that Y takes on values much larger 
than a. 

To fix this one might try 



C (5) = C (l) -a, i = 1,2,3 or 4, (2.20) 

since the mean E(C^) is easily calculated if E(C^) is 
known for i = 1,2,3, or 4, and E(a) is known. But E(a) is 
l/(l-p) for any M/G/l queue (Cohen, 1969) ; approximations for 
the GI/M/1 case are discussed in the Appendix. 

There is an additional problem of dimensional stability 
involved in using , as with C* 1 * , , in that 

E(C^) depends on both y and p, while E (a) depends only on 
p. Thus control is not uniform across the whole range of param- 
eter values. 
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To avoid this one use 



C (6) = C (l) - a/y. (2.21) 

This, with i = 2, was found, in simulation studies, to be 

the most successful control variable in that it obtained a variance 

reduction which was uniform in y and p (or y and X) and its 

mean value is fairly simple to compute. These simulation results 

are discussed in the next subsection. 

Note that multiple control variables using any of the above 

controls can be used. In particular, one need not take the differ- 

( 2 ) 

ence of say, C and a/y, but may use a multiple control. 

However, the fact that two regression coefficients, say 6^"*^ 

( 2 ) 

and 3 t must be estimated from the simulation data makes the 
possible gain in variance reduction of dubious value. 

Note also that since all the control comes from within 
cycles, there is no reason that the confidence interval estimation 
techniques referred to in the introduction would not go through 
for the variance-reduced estimates. This has been verified in 
simulation runs. 

2 . 3 Internal Control: Simulation Results 

It is not possible to verify analytically what variance 
reduction will be obtained via the several internal controls 
listed in the previous section, or to get an idea of the magnitude 
of the effect. Even for something as simple as it is diffi- 
cult to compute analytically the correlation between and 

Y-a for the M/M/1 queue, and this is what is required in the 
equation (2.5) to find the variance reduction. 
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Thus, we resorted to simulations to verify the amount of 
variance reduction obtained and the relative effectiveness of the 
various controls. In the final simulations all runs were performed 
on IBM System 360/67 computer using the LLRANDOM package (Learmonth 
and Lewis (1973)) which generates random numbers according to the 
scheme given by Lewis, Goodman, and Miller (1969) and exponentially 
distributed random numbers using the Marsaglia "rectangle-wedge- 
tail" method. Tests of the random number generator are given in 
Learmonth and Lewis (1974) . 

Of the extensive simulation checks performed, we give here 
only a summary of the conclusions and one detailed tabulation and 
one short tabulation in the case of the most suitable control. 

(1) The controls and do much better gener- 
ally than C^, with little improvement over obtained 

by use of and C^. We say generally because results 

vary with X and y and their ratio p. 

(2) Subtracting the number of customers served in a busy 
period generally improves the variance reduction. By making 
it dimensionally stable as in (2.21) with i = 2 we obtain a 
"variance reduction" measured in terms of ratios of standard 
deviations, of approximately 70%, uniformly over X and y. 

This is roughly equivalent to halving the number of cycles (b.p.'s) 
that one must simulate; (0.7) 2 ~.5. Much better reductions can 
be obtained for smaller p (i.e. p = 0.25) by specially designed 
controls; the point is that C ^ ^ using C ^ ^ works even out 
at p=0.99, where variance reduction is extremely important. 
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Table 1 shows results obtained by simulating an M/M/1 queue 
out to n = 2000 cycles and replicating the simulation 250 times to 
estimate the variance of the estimates W(n), W CT (n), W CB (n), 
where we drop the n for convenience. Here, we have specifi- 
cally that 

1 
n 

W (n) = , 

k X ta k+ 6 B (C (6) -E{C <6) })] 
k=l 

where 

C (6) = C (2) - a/y. (2.22) 

The estimated precision (standard deviations) of the estimates of 
E (W) are given in brackets under the estimates. 

The results in Table 1 are for p = 0.5 and three values of 
y; the results are not very different at different values of p. The 
case p = 0.99 is given in Table 2. The second, third and fourth 
columns in the Tables give correlations between the control and 
Y-a etc., from which the theoretical variance reduction can be 
computed. They are very close to the values given in the next 
to last column, from which we deduce that estimating $ T and 
affects the variance reduction only slightly. There is negligible 
effect of different values of y for fixed p = 0.5 and fixed 
P = 0.99. 

Note that for the results for p = 0.99 given in Table 2, the 
variance reduction is 73% (about the same as for p = 0.5). For the 
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case where the control is on the top, i.e. for Y, the variance 
reduction is not quite as good for control of a. Note too that 
the estimated values appear in some cases to be at least three or 
four standard deviations from the true mean. This is because the 
estimates W, W T and W g can be seen from the 100 replications 
to be non-normal. In other words, for high p(0.99), the simula- 
tion needs to be taken out further than 2000 cycles. 
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3. Internal Stratified Sampling 



Another technique for variance reduction which can be 
potentially used with the regenerative method is stratified sampling. 
For a brief description see Kleijnen (1974) , p. 110. In essence 
this uses analytical information in the following way. 

If we can stratify or partition a random variable X by 
its values (or those of a concomitant variable) into K strata 
labeled k=l,2, ,K, we can write the mean of X and the sample 

mean X as, respectively. 



K 

y = E ( X) = l P, E ( X | Xe strata k) ; 
k=l 



X 



1 n 
1 l 



X. 

n i=l 1 




= l 

str 1 




n 




+ l 

str K 




n 




/ 



(3.1) 



(3.2) 

(3.3) 



where n, = number of X. 's observed in strata, and p, is the 

probability of being in strata. 

Now, if the p^ ' s are known and we substitute them in (3.3) 

for n^/n, we 9 et X , a stratified estimate. It will be biased, 

since the divisions of the sums in the populations are random; if 

the numbers observed in each population are controlled and taken 

to be np^,np 2 » etc., we have what is called a proportionally 

sampled estimate X with 

ps 



var(X ) 
ps 



K 

var(X) - l P (y -y) 2 /n, 
k=l K k 



(3.4) 
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where y^ = E{X | X e strata k}. Thus, because of the use of prior 
analytic information we always get variance reduction. 

The variance of X gt is not analytically tractable, but 
early studies reported in Kleijnen indicate it is close to (3.4) 
if n and all the n j c ' s are sufficiently large. Our simulation 
studies with stratifying a have confirmed this; because of the ease 
of computation of P{a=l} ,P{a=2} ,P{a>3} , as in section 2.3, it is 
natural to stratify a. Considerable variance reduction is obtained, 
especially for smaller values of p, the traffic intensity. 

Unfortunately, when the quantity a is stratified in the 
bottom (denominator) of the estimator W(n) very little overall 
variance reduction is obtained unless p is small. We do not 
understand this lack of variance reduction but apparently the 
correlation between the stratified version of a(n) and Y (n) 
is reduced. Analytical studies of this effect are very difficult. 

It is also possible to use a stratified estimate of a 

in W(n) and to control the top, Y (n) , with say C v ^ using 
( 2 ) 

C . This works well for small p, but increases the variance 

as p approaches 1. Again it is difficult to understand this 
effect. 
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4 . Conclusions and Summary 



We have been able to obtain a worthwhile variance reduction 
using internal control variables, for the regenerative estimate of 
the limiting value of the mean waiting time in an M/M/1 queue. 

This reduction is obtained uniformly over all parameter values. 

It is fairly certain that the technique will work well with any 
GI/G/1 queue or other regenerative stochastic processes or systems. 
Internal stratified sampling schemes, however, did not work nearly 
as well. 

The techniques can be extended to other stable stochastic 
systems, such as the Markov chains considered in Crane and Iglehart 
(1974b) . In that case the computation of the mean values of the 
controls is simpler because of the structure of the Markov chain. 

The main problem in applying the internal control variance 
reduction technique seems to lie in the fact that the estimator 
proposed by Crane and Iglehart (1974a) involves a ratio of two 
random variables, and these are difficult to work with in general. 

An alternative which will be considered later is to use the 
existence of regeneration points more specifically to obtain vari- 
ance reduction with the classical estimator W given at (2.1). 

m 

One advantage which the regenerative estimator W(n) has over 
W m is the ease of obtaining confidence interval estimates or esti- 
mates of the precision of W(n) and W(n). This is not a draw- 
back if the simulation is large and more than one (say ten or 
twenty) realizations of the queue are obtained. 

To fix ideas note that we can write W as 

m 
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m 



J Y + Y ' 

l k=1 k + *N(m)+l 




where N (m) is the number of completed busy periods in the queue 
in [0,m], Y^ as before is the sum of the waiting times in the 



incomplete cycle. 

Now it is possible to apply internal controls to each Y^ 

in the sum. Problems arise in estimating the coefficients 8 in 

the control because they involve a random sum of random variables. 

But it is much easier to find a control C for Y. rather than 

k 

the difference Y^-a^, and also it is still possible to apply 
external controls as well as internal controls. 

These ideas will be followed up in a later paper. 



k 



th 



cycle, and Y N( m )+i t * le sum °f the waiting times in the last, 
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APPENDIX 



To implement either the internal control or stratified 
sampling techniques certain theoretical parameters associated 
with the GI/G/1 queue are required. In this appendix we shall 
indicate the values of these parameters in so far as they can be 
calculated. These values are either well-known or easily calcu- 
lated. For a reference to the known formulas see Cohen (1969) . 

We begin with E{a^}, the expected number of customers 
served in a busy period. For the general GI/G/1 queue recall 

that we let X = v - u and S = X, + . . . + X , for n > 1, 
n n-1 n n 1 n 

with Sq = 0 . Then a^ = inf{n>0: S n <0}. The general expression 
for E{a^} is given by 

OO 

Efa^} = exp{ £ n -1 P[S n >0]}, 

n=l 

an impossible expression to evaluate in general. Another useful 
expression for E{a^} is 

E{a } = 1/P{W=0}, (A . 1 ) 

where W is the stationary waiting time. In the special case of 
M/G/l, however, we have 

Eto^} = (1+p)" 1 . 

Now for the queue GI/M/1 we can use (A.l) and the stationary 
distribution of the embedded Markov chain to conclude that 

E {a 1 > = ( 1-6 ) -1 , 
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where 6 is the root inside the unit circle of 



z - U{y (1-z) } = 0 

" su l 

with U(s) = E{e }, Re s > 0 , and u^ is an exponential (y) 

r.v. It is easy to check that 6 = p for M/M/1 queues. Daley (1975) 

has recently proposed the approximation to 6 given by 

S = a^ (1-p) 2 + 2 (1-b -1 ) p + (2b -1 -l) p 2 , 

where a 1 = P{u 1 =0} / E{u 1 > = 1, and b = E{u^}. This approximation 

gives good results in a number of examples calculated by Daley (1975) 
and may be useful for the purposes of this paper. 

Next we turn to the computation of P{a^=l) and P{a^=2}. 

For the GI/G/1 case we have 

P{a 1 =l> = PtSjSO} 

and 

P{a 1 =2 } = P{S 1 >0,S 2 <0}, 

both of which can be worked out with a little effort. For the M/M/1 
queue 

P{ a;L =l} = (1+p)" 1 

and 

P{a i=2 } = p(l+p)" 3 . 

For the M/G/l queue 

P{ a;L =l} = V ( A) 

-Av n 

where V ( A) = E (e ) , and for the GI/M/1 queue 
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P{a 1 =l} = 1 - U(y) , 

where U(s) is given above. For the M/E^/l queue and the E^/M/l 
queue the value P{a^=2} can be calculated with some effort. As 
these expressions are cumbersome they shall be omitted. 

Next we turn to various partial expectations which are 
needed for internal control 



E{S l +S 2' a l 



a 1 >2} = E{S 1 ,S 1 >0} + E{S2fS 1 >0} 



and 



E{a 1 ,a 1 >2} = E{a 1 } - P{a 1 =l} . 



In the special case of the M/M/1 queue 



E{S 1' S 1 >0} = * 




and 



E{a 1 ,a 1 >2} = 2p/ { ( 1-p) ( 1+p) } . 
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